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("-^ , By generalizing the usual current density to a matrix with respect to spin variables, a general 

fvj ■ equation of continuity satisfied by the density matrix and current density matrix has been derived. 

This equation holds in arbitrary spin-orbit coupling systems as long as its Hamiltonian can be 
expressed in terms of a power series in momentum. Thereby, the expressions of the current density 
matrix and a torque density matrix are obtained. The current density matrix completely describes 
both the usual current and spin current as well; while the torque density matrix describes the spin 
OO ' precession caused by a total effective magnetic field, which may include a realistic and an effective 

O^ . one due to the spin-orbit coupling. In contrast to the conventional definition of spin current, this 

expression contains an additional term if the Hamiltonian includes nonlinear spin-orbit couplings. 
Moreover, if the degree of the full Hamiltonian > 3, then the particle current must also be modified 
in order to satisfy the local conservation law of number. 
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I. INTRODUCTION 



Pit Pll 



Recently, the spin transport, precession, and accumulation generated by spin-orbit coupling and the resulting 

intrinsic spin-Hall effecliii^ have attracted considerable interest by physicists^'^'^'^'^'^'^ i^^'^-'^d^i"'^^ because it may provide 

i-pj ' a promising method to manipulate the electronic spin without application of magnetic fiel d-'^'^d^'^^'^^ in spintronics i^^d^ 

r^ I However, some conceptual ambiguities still remain in how to properly describe the spin transport and precession in 

O ' spin-orbit coupling systems. In order to describe these phenomena, a conventional spin current operator is intuitively 

C^, defined as Jl = {l/2){vi,Sj}, with Vi — dH/dpi a velocity operator and Sj is a spin component (actually, as will be 

pointed out in Sec. IV, this expression should also contain a factor of density operator (S(x — x) in order to have the 

00 ] dimension of a current density). However, it is still very questionable whether this quantity is meaningfully observable. 

t> ' This problem has caused a lot of confusion and debates in recent literatur o^°'^^i^^'^'^i^'^'^^ and has even become an 

^^ ^ obstacle to the further progress of this filed. From the experimental point of view, the most interesting observable is 

vN , the spin densityr^ii^iii which can be given by tr{sjp), with 

o 
\o 
o 

the density matrix; i.e., the density matrix p can completely describe the number density as well as its spin polarization. 

d ' Theoretically, a fundamental equation describing a transport phenomenon of any additive conserved quantity, say Q, is 
the equation of continuity dpq/dt + V • jg = 0, which connects the change of its density pq to a corresponding current 
I ' Jq; if this additive quantity is not conserved, we must introduce a source (or drain) term Sq in this equation such that 
' O \ dpq/dt + W -'^Q = Sq. As a classical picture of motion, the spin transport in spin-orbit coupling systems is very similar 
to the latter case, since the spin is also an additive and local quantity (the locality means that the spin commutes with 
the position); i.e., the spin density can be defined as a function of position. Therefore, a straightforward approach to 
introduce the spin current is to establish a continuitylike equation satisfied by the density matrix p. This equation 
should be a general counterpart of the usual continuity equation of number dp/dt + V ■ J = 0, and can be called 
a general continuity equation. We expect that a corresponding current density matrix J can be naturally obtained 

$H from it. However, as the spin is not conserved, so the time rate of change of p must contain two distinct parts, i.e., 

dp/dt — —V • J + T. The first part is a divergence —V • J, which can be ascribed to a current density matrix J; 
while the second one, denoted by T, does not possess this character and can be called a torque density matrix. If 
we integrate both sides of this equation over a volume 51, then we have (d/dt) J^ pdV = — §q^ dS ■ J -I- /^j TdV . The 
surface integral ^„„ dS ■ J can be interpreted as a current flowing through the boundary surface 9f2, while the volume 
integral J„ TdV represents the precession rate of the spin polarization within the volume fl. As will be demonstrated 
in Sec. II, this division can be uniquely determined by the equation of motion dp/dt — [_ff(p,x),p] = £[p(x, x';t)] 
and the condition that T can no longer be expressed as a divergence, where C[p\ = >C(V, V')[/5] is a functional operator 
acting on p and V — Vx and V' == Vx' are the differentiation operators with respect to x and x', respectively. The 



essence of this problem is to cast the C[p\ into the form 

L{V, V')[/5] = (V + V) • J[p\ + r(i(V - V))[p], 

where J7'(V, V')[p] and T(i(V — V'))[p] are also functional operators. It will be proven that J and T can be uniquely 
determined by the equation of motion as long as C can be written as a power series in V and V'. So, if we define 

J(x) = - lim j'(V,V')[/5(x,x')], 
f(x)= limT(i(V-V'))[/5(x,x')], 

x'— »x Z 

then the continuity-like equation can be obtained. It is important to note that if we transform the /5(x, x') to the 
Wigner function /5(p,R), where R = -^(x + x') and p is the canonical momentum corresponding to x — x', then 

the expression for T has the form r(R) = /d'^pT(ip)[y6(p, R)]; it no longer contains the spatial derivative Vr and, 
therefore, can not be written as a divergence. 

From this general continuity equation, we can deduce that the trJ and tr(sJ) represent the current density of 
particle J and the current density of spin J^, respectively, just as the trp and tr(s/5) represent the number density p 
and spin density Ps(x), respectively, i.e., the matrix J can provide a complete description of the particle transport as 
well as the spin transport. In order to confirm that this physical interpretation is self-consistent, we must prove that 
(1) trT = 0, which is a necessary and sufficient condition for the local conservation of number and (2) J and T must 
be gauge invariant in order to ensure that they can be interpreted as physically observable. In this paper, a general 
equation of continuity satisfied by p, J and T is derived, which holds in arbitrary spin-orbit coupling systems as long 
as its Hamiltonian can be written as a polynomial in momentum. Consequently, the expressions of the J and T are 
presented. These expressions are gauge invariant and independent of the magnitude of spin and the system dimension 
and satisfy trT = 0. In contrast to the conventional definition, this formulism contains the following corrections: First, 
if the spin-orbit coupling Hamiltonian includes second- or higher order power of momentum, e.g., in the Luttinger 
model,— 1^ the two-dimensional (2D) cubic Rashba model,— or the Dresselhaus spin-orbit coupling , ^"^'^^i^^ then we 
should add a new additional term to the expression for J, which is purely originated from the nonlinear spin-orbit 
coupling and contributes only to the spin current. Second, if the degree of the full Hamiltonian degi?(p) > 3, e.g., 
in the cubic Rashba model or Dresselhaus spin-orbit term, then the particle current density J as well as the spin 
current density are also different from the conventional definition; otherwise, it will violate the local conservation law 
of number. Moreover, this equation of continuity also holds in the presence of arbitrary electromagnetic field and 
Coulomb interaction or nonmagnetic impurity scattering. Consequently, from this continuity-like equation we can 
also deduce some exact identities which are very similar to the Ward-Takahashi identity, although the spin is not a 
conserved quantity. 

This paper is organized as follows. In Sec. II the general equation of continuity in a spin-orbit coupling system 
without Coulomb interaction and impurity scattering is derived. Then, we will prove that this equation also holds 
in the presence of the Coulomb interaction and nonmagnetic impurities in Sec. III. Section IV is devoted to a detail 
discussion and comparison of these expressions with the conventional intuitive definition. A short summary of the 
conclusions will be given in Sec. V. 

II. GENERAL EQUATION OF CONTINUITY 

Firstly, we consider a general spin-orbit coupling system without any two-body interaction or impurity scattering, 
but may subject to an external electromagnetic field described by {</>, A); the case in the presence of the Coulomb 
interaction or impurity scattering will be discussed in Sec. III. The corresponding one-body Hamiltonian has the form 

H = -— (P--A) +ecj) + Ko = hn + ecj) + Ko, (1) 



2m V c 

where h^ = (1/2to)(P — -A)^ is the kinetic energy and fiso is a general spin-orbit coupling Hamiltonian and we 
suppose it can be expanded as a power series 

^so = -ff<°^ + H\^H^ + H\fp,pj + Hlf^P;PjPk + ■■■ ■ (2) 



Here and hereafter, a summation over repeated sufRx is implied, where pi = —ihdi — {e/c)Ai = —ihdi {i = x,y, z), 
di are the covariant derivatives, and H^',,^ {n = 0,1,2,---) are the operators acting on spin variable (or 2x2 
matrixes if the spin is 1/2); meanwhile, i?^"-* is a scalar representing Zeeman energy due to magnetic field, H^ are 
the components of a vector, H^ ..^ are the components of a tensor of rank nth and are supposed to be symmetric 

^ (n) 

under the permutations of its suffixes ii,Z2, . . . , ««. Furthermore, we assume that all H^ ...^ are independent of 
position except ij'"'. The corresponding many-body Hamiltonian is H = /^ct(1)^ctct''0ct'(1)i where a and a' are spin 
variables. 

In order to derive the general equation of continuity, we consider the equation of motion of Green's function, for this 
approach is apparently gauge invariant. The density matrix p can be given by pap^xt) — — (=^[?/'a(xi)'!/'/3(x<^)]) = 
—ihGapi^t, xt"*"), where £^ is the time-ordering operator and t+ = t + 0+. From the equation of motion wc have^i^ 



ih^-ecj,{l))-(-th^~e^{l') 



G(l,l')- ho + ko,Gil,l') 



-0, 



(3) 



at ^' 7 V dt 

where G denotes the matrix Gap, and (1) = (xi), (1') = (x'i'). By taking the limit as 1' -^ 1+, the first term becomes 

dpi^t) 



lim 
1'— 1+ 



ih— ~ e0(l) + ih— + e0(l') 



dt ^' ' at 

The commutator [ho, G] in the second term of Eq. ([3]) can be written as 



G(l,l')-- 



at 



(4) 



ho,G 



-^{P--A,G(1,1')} 



c 2m I c 



where (P-(e/c)A)G'(l,l') = [-inV-(e/c)A(l)]G(l, 1'), G(l, l')(P-(e/c)A) = [iW-(e/c)A(l')]G(l, 1'), [X,Y] = 
'E^iX^,Yi\, and {±,Z} = J2^{^^:Z}. Bccausc 



hm |P - -A, G(l, l')| = ^ |P - -A, p{l, l')| 

1'— ►!+ L C J n L C J l' = l 



SO, if we define a current density matrix 

Jo(l)^-^{P--A, /5(1,1')} 

Zm L c J i'=i 



then we have 



lim 



fto,G(l,l') =V-Jo(xi). 



(5) 



(6) 



Now, we consider the spin-orbit coupling term [hso, G(l, 1')]. According to Eq. 

1 r. 



lim 



7J("),G 



ij("), 



i'=i 



the first term is simply 



(7) 



In order to treat the other terms of [/iso,G(l, 1')], we note that H^',^ (n > 1) commute with pi because they are 
independent of the position, and use the following operator relation: if [A, B] = 0, then [AB, C] — {1/2){A, [B, C]} + 
(l/2)[i, {B, C}]. By substituting Eq. ^ into [hso, G(l, 1')], we obtain 



"■so, ^ 



I E {{h^:U . [pn ■■■p..,g]} + [ij("'.,„ , {p, . . . p.„ , g}] ) 



(8) 



where 



Ki ■ ■ ■ P»,. , G 



[Pn■■■p^„-pl■■■pl]G{l,l'), 



(9) 



{P^, ■ ■ -P*,., g} = {p^, ■ ■ -p,^ -t-4 • ■ -p^J G(l, 1'), 



(10) 



and p' = ihd'j — {e/c)Aj{l') = ihd'j, which commute with pi because pi act on 1 = (x, t) while p' act on 1' — (x',t'). 
It is easy to verify that Eq. ([7]) is also a special case of identity ([5]). The limit of the first term on the right-hand side 
of Eq.® can be entirely written as a divergence form. For n = 1, we have 



^lim^ i{i7«, [p,, G]} - V • i {h(i),p(1, 1')}^,^^ 



where H*^^^ — [Hx , Hy , Hz )■ For n = m + 1 > 2, we define an auxiliary operator 



(11) 



Av- ,„(A:,pO = Pii ■■■Pz^+Pzi-- ■Pi^.-iPL + • • • +Pii • • -pL = JZ^^i ■ ■ -P^^-iPL-i+i 



■■■Pi. 



1=0 



An) 



Owing to the symmetry of -ff , " ...j and commutability between pi and p', we havi 

1 



,32 



So, if we define 



{^X^l..' [pnP^. ■ ■ -P.™.., G(i, 1')] } = iP.-pd^ {h^Z^d^I^Jh, 1')} 



(12) 



(13) 



(14) 



where we also introduce a D^^' = 1 in order to include the term of Eq. (jlip in this unified formula, then we have 

1 



1,1™+ E 2 {^^--"' [P^^ ■ ■ ■P^r.,G{i, 1')] } - dA) 



= V • Jso. 



(15) 



The limit of the second term i A- " .i„ : "^ Pn " ' ' Pi„ G M on the righ-hand side of Eq. ([5]) cannot be entirely written as 
a divergence. However, we can prove that it can be written as a divergence and a remaining term, while the remaining 
term no longer contains spatial derivatives Vr if it is represented by the Wigner distribution function. In order to 
prove this conclusion, we take first the case when n = 1, it evidently cannot be written as a divergence as 1' — » 1^, 
because {pi, G} = —ih(di — 9j')G(l, 1') and only the limit of {di + 9,-) can be written as a divergence. Then, we take 
n = m + 1 > 2 and prove that the limits of 



H^ 



:^i...„{a.p.---a„,g} 



and 



[p^lP^2■■■P^,.,G]H^^'^l■■■^„ 



_2fj(n) Pn +PnP^2+Pi2 P2^_±p1^a 
iii2---j„ 2 2 '" 2 



r,P^l + P'h P^2 + pL P«" + pL n ff (") 

^ ' '^^M«2-i„ 



2 2 2 

can be written as a divergence form. To this end, we again introduce the following auxiliary operators: 



?(") 



K'-^.iP^^P'i) 



Pll + Pn P^2 + Pt2 Pir. + pi 



Gi^-l^PriPi) = Ph ■■■Pi^ +Pil ■ ■■Plrr^-lR 
- pi ■■■P'lrr. ~Pi ■■■PL-1^ 

and define another term of the current density matrix J^^ as 

(^o).-Ei[^"^l'C^"),„/5(l,l' 

m>l 

and a torque density matrix as 

^[^]-Ei[^::l-^„'^l:t...„p(i'n 





2 


; 






(1) 




+ P^ 




-1) 


(1) 




-Pi 




-1) 



(16) 



(17) 



(18) 



(19) 



It is easy to verify that the terms of n = 0, 1 can also be incorporated into the righ-hand side of definition p9|l if we 
define M^^ = 1. Then, we can prove that (see the Appendix A) 



i/™+ Yl 2 [^lt-»„' {p^iP'2 • •■?.„, g} 



V-3i,-f[p], 



(20) 



where J^^ is a new term of the current density matrix, which is originated from the nonlinear spin-orbit coupling, 
because it will vanish if fl^j j ...j are independent of spin (i.e., it is only a c-number rather than a spin matrix). More 
importantly, we should note that the T[p\ cannot be rewritten as a divergence because it is given by a polynomial 
of di — d[. This property will be more evident if we transform the p(l, 1') to a gauge-invariant Wigner distribution 
function^^ by 



P' 



^„Mp,£;R,t) = |d3rdrei^[(-~+^*(^*)l-[P+f^(^*)l-'->(^t(R_ |^,„ ^)^^(R+ £^,+ ^))^ 



in which the macroscopic variables R = (x + x')/2 and microscopic variables r = x — x' have been used, and then we 
have 



H, 



(») 






(27rfi)4 






,Pii •••?»„, /5(p,e;R.;0 



(21) 



It does not involve the macroscopic spatial derivative Vr, and therefore can not be written as a divergence Vr • J. 
Comparing this expression with Eq. ([7]), it is obvious that the H^ ,,^ pi^ ■ ■ -pi^ plays a role which is very similar to 
Zeeman energy H^'^' . The only difference is that the former depends on momentum, because it is originated from 

^ (n) 

the spin-orbit coupling, while the latter depends on position. So, the H^_^,,,^^pi-^ ■ ■ -pi^ can be regarded as the energy 
contributed by an effective magnetic filed depending on momentum. Take the limit as 1' — > 1"*" and substituting Eq. 
(HI), (O, (HH) and (pOl) into Eq. ^, we finally obtain the following generalized equation of continuity: 



dp 



V- J 



Jso =r[p]. 



(22) 



Evidently, the form of this equation is independent of the system dimension or the magnitude of spin. The only 
differences are the number of components of the current density and the dimension of the matrices, if the system 
dimension or the spin magnitude are different. 

From the definitions ^, (fT4|) .(|18 p . and ([T9|) . we can deduce that J and T are gauge invariant because all the 
derivatives in their definitions are covariant.^'^-'^^ Moreover, from the definitions (|18p and (|19p . we get trJ^^ = and 
tr(T[/5]) = 0, owing to their commutator form. The identity tr(T[/5]) = ensures the local conservation law of number 
(or charge), while trJ^Q = means that the J^q does not contribute to the particle or charge current, but it may 
contribute to the spin current. Because tryo(xt) is the number density /9(xi), so, trJ(xt) = tr(Jo -I- Jso) is the current 
density of the particle J(xt), while the charge current is Je(x<) = eJ(xt), here the identity trJ^^ = has been used. 
By letting Jq = trJo and Jso = trJgo and taking the trace over spin of both sides of the Eq. (j^^ . we obtain the 
conservation law 



| + v,p. 



= 0. 



(23) 



Similarly, because 



Ps(xi) = tr[p(xt)s] = tr [sp(xi)] = -tr {p, s} 
is the expectation value of the spin density, therefore, it is natural to define the spin current density as 



J'*(xt)=tr J(xi)s =tr sJ(xt) =-tr|j,s| 



and the spin torque density as 



T^(xi) EE tr f(xt)s = tr sf(xi) 



-tr 



{t.} 



(24) 



(25) 



(26) 



From Eq. (|22p . we obtain the following equation of continuity satisfied by the spin density, spin current density and 
spin torque density: 

^ + V ■ r = T^ (27) 

This equation and Eqs. (flQ)) and pT|) reveal that the matrix T completely describes the spin precession caused by 
a total effective magnetic field, which may contain both the Zeeman term and the effective magnetic field originated 
from the spin-obit coupling. 

A potential question is that whether the torque density matrix T can also be rewritten as a divergence by solving 
the equation T = — V • JTrr^ although the expression for T itself no longer explicitly contains spatial derivatives Vr. 
This is not practicable, because the condition T — —V • 3t is not sufficient to uniquely determine the expression for 
3t except in one-dimensional systems, since the 3t will have two or three unknown components in a 2D or three- 
dimensional systems, but there is only one equation, which is evidently not sufficient. ^^ The physical meaning of this 
conclusion is that the precession rate J^ TdV in a volume O cannot be determined only by the information on the 
boundary surface dfl, i.e., it cannot be expressed by a surface integral such as a current. In contrast, the effect of 
the term V • Jg^ can be expressed as a surface integral Jg^ J^q • dS if we consider {d/dt) J^ pd'^R. Furthermore, from 
definitions ([TT]) and p^ . we can find that there is a common factor pi—p'^ = —ikdn- in the expression for J^^, which 
implies that the term Jg^^ J^^ • dS may still be nonvanishing when even there are no any electrons inside the volume 
ri, i.e., /5(R) = for all R G il, because the derivatives dp/dRi on the boundary may include a nonvanishing spin 
polarized part. So, it may give a nonvanishing J^^, and the {d/dt) J^ pd^Ti ~ — L,, J^^ • dS may also be nonvanishing; 
it turns out that the time evolution rate of the total spin inside the volume is nonvanishing even when there are no 
electrons. Therefore, the /g^ Jso ' ^S can only be interpreted as another kind of current fiowing from the outside of 
the volume f2, since in this case there are no spins (or electrons) inside the volume; so, it is impossible to be explained 
as the spin precession inside the volume. 

III. THE INFLUENCE OF THE COULOMB INTERACTION AND THE SCATTERING BY 

NONMAGNETIC IMPURITIES 

Now, we consider the infiuence of the Coulomb interaction, which can be written as 

H, = ^ J dldl'i,lil)^jliV)v{l - l')ij.<{V)M^\ (28) 

where v{l — 1') is independent of spin. So, the total Hamiltonian is H + He, and Eq. fS]) become o*^^'^^ 

[Go iG - GGo i]„^(l, I') = -ihj[v{l - 2) - v{l' - 2)]G„^,^Mir , 1'2+), (29) 

where [Go"^G](l, I') = ih{d/dt)G{l,V) - h{l)G{\,l') and [GGo"^](l, I') = -ih{d/dt')G{l,V) ~ G{\,l')h{V)- 
Ga-y.-y'/3'(12, 1'2') = [l/(i?i)^](T[f/;a(l)?/'^(2)V'!,/(2')V'L(l')]) is the two-particle Green's function, and a bar over vari- 
able (2) indicates that it is the integral variable, 2~ = (^2,^2 — 0+). As V -^ I+, the left-hand side of Eq. ((29l) is 
T[p] — (dp/dt) — V • J, while the right-hand side will vanish because 

[v{l-2)-vil'~2)]r=i = 0. 

Therefore, we again obtain the Eq. ((22|) . 

We take into account the scattering of some random nonmagnetic impurities and suppose that the impurity potential 
has the form [/(x; {ri}) = J^i'^i^ ~ ^i)^ with r^ is the position of ith impurity and here U{{ri}) is a function of 
the set {r.;}. The additional term of the second-quantized Hamiltonian Himp({ri}) = J dlU{l;{ri})'il;i{l)^p{l). 
Now, the interesting quantities become some quantities averaged over these random impurity positions, e.g., p(l) = 

/ • • ■ JYl^{d'^ri/n.)p{l; {vi}), and J(I) = / • • • JYl^{d^ri/V,)3{l; {vi}), etc., where ^ is the volume of the system. So, 
the Eq. ^ becomes 



ih— - e(j){l) + ih— + e(biV) 



GdrJ)- h + U{{r,}),G{{r,}) =0. (30) 



Because U is independent of spin, we have 



JunjU{{r,}), G(l, l'; {rj)] = ^^JU{1, {r,}) - U{l', {r J)]G(1, l'; {r J) = 0. 



Therefore, we get 



V-J(xt;{ra) = f[p({ra)]. (31) 



dt 
By taking the ensemble average withe respect to {r^} of both sides of this equation, we have 

i i 

interchange the ordering of the integrations and differentiations, we obtain 

g + V.J=T[;5]. (32) 

Similarly, it is can also be proved that this equation still holds if the Coulomb interaction and the impurity scattering 
exist simultaneously. Furthermore, if we take the functional differentiations with respect to the external fields on both 
sides of this equation, we will obtain some exact identities connecting the correlation functions of the spin current, 
spin torque, and the usual particle current, which will be very similar to the usual Ward-Takahashi identity even 
though the spin is not a conserved quantity. So, if we use some approximation methods to calculate the response 
of the spin current or spin torque to external perturbations, then the self-consistency must be taken into account 
in order to preserve these exact identities, just as in the problem of the electron transport. '^^"^^ In addition, these 
expressions can also be applied to the finite temperature systems, only if all averages are defined with respect to a 
grand canonical ensemble, in which the spin-orbit coupling terms will appear as a thermodynamic weighting factor 
exp(— Hso/fcT). Because the spin-orbit coupling coefficients are often very small, so Green's functions and the related 
observables may appear as spin polarized only if the temperature is low enough such that Hgo/kT can not be ignored. 
Otherwise, this weighting factor will be almost spin isotropic, so the phenomena due to the spin polarization will be 
unobservable. 

IV. DISCUSSIONS AND COMPARISONS 

The expressions derived in Sec. II, which are expressed in terms of density matrix, can also be given by the 
corresponding one-body operators such that (/) = tr/dx(ix'/(x',x)p(x, x') = Tr(//5), where/(x',x) = (x'|/(p,x)|x) 
represents a one-body operator; "Tr" denotes the trace over both the spin and coordinate variables (while the "tr" is 
the trace over spin). This form is more convenient to compare with the conventional definition, since the latter was 
usually given by a one-body operator. To this end, we need only to consider a special case /5(x, x') = ^(x)?/'*(x') and 
obtain the following relation: 

[■P^^ ■ ■ ■Ptrr.P'j^ ■ ■ •Pj„/5(x,x')]x'=x = {lp\Pn ■ ■ -PjJi^ - x)P»l • • •-P»„JV') 

= ( —ihdj^ ) ■ ' ' ( ~*^^i.i ) i'i^) ( —ifidii ) ' ■ ' ( ^'''^'^im ) V'(x) 7 (33) 

where x is the position operator and x is the eigenvalue of x [note that 5(x — x) — |x)(x| is the probability density 
operator]. Here, the spin suffixes are ignored for simplicity. From this relation, we have the following corresponding 
rule: 

Pn ■ ■ ■ PtmP'ji ■ ■ ■ P'j„ ^=^ Pji ■ ■ ■ P]J{^ - x)Ki ■■■Pi^- (34) 

According to this rule, we can obtain the corresponding one-body operators of Jq, J^, Jso, Jgoi (JD'j ^^'^ '^^j 
which will be denoted by jo, Jq: jso, Jso, ^^^ ^^' respectively. From definitions ^ and Jq = trJg, we have 

Jo = :5^{p,<5(i-x)} (35) 

2m 



and its mean value (V'lJolV'} = {ih/2m) (V'i/;)*V' ~ V"*^^ • Its contribution to spin current density is 



In order to express the one-body operators jso and jg^, we define 

d-"..j^(p,x) = (5(x-x)pij •• -pi^ +pi„^(5(x-x)pij •■ •_Pi„_i H \-pt^ ■■ •pi„/(x-x). 

This is the one-body operator corresponding to definition (fT2| . According to definition (|14p . Jg 
JgQ = tr(sJso), we obtain that the term of particle current density jso can be given by 

while its corresponding spin current density can be given by 



(36) 



(37) 



trJso, and 



(38) 



(39) 



Similarly, let 



and 



Am) 



2 ' 1 2 



-(™) 



-,,5(x-x) 



(l)l 



From definitions ^, dH]), ([251) and JM]) we have 



I [Ai: 'i2---i„ 



(40) 



(41) 



(42) 



i" = 



^i['' 



H-(") 



n*2- 



(43) 



In order to compare these expressions with the conventional definition, we must first clarify a misleading statement 
in the conventional definition. In recent literature, the spin current was usually defined as J^ = 2 {'^, s}. However, what 
we are actually concerned is not the mean value of this operator ('0|J^|'0), but the mean value of some kind of density 
of this operator i[V'^(x)(J**^(x)) -I- (J**-0(x))'''f/;(x)] — (V'l^jJ^, (5(x — x)}|?/;). Therefore, more precisely, the quantity 
we actually defined is |{{v, s}, (5(x — x)} rather than J** itself and should be called a spin current density. Moreover, 
the spin current, by analogy with any other kind of current, should be defined as an integral /**[A] = Jj^j^dAi of the 
spin current density jf over a surface A, which is also not the |{v, s} itself. Now, a puzzle that confronted us is that 
there exist some different Hermitian combinations of the products of v, ^(x — x) and s because they do not commute 
in spin-orbit coupling systems. They all have the same dimension and the same classical analogy (or classical limit), it 
is very difficult to determine which one should be interpreted as the proper spin current density just by this intuitive 
method. For example, we can also define |{v, {(5(x — x), s}} or |{{v, S{x — x)}, s} as a spin current density (only the 
first one and the third one are equal owing to the fact that s and x commute). From the classical picture of motion, the 
physical meaning of the second one can be interpreted as the velocity multiplied by the spin density; while third one 
is the particle current density multiplied by spin. Moreover, if the velocity includes higher power of momentum, i.e., 
it is a nonlinear function of momentum, then there will be even more probable expressions of the operator due to the 
noncommutative property of the density (or position) and momentum. According to the results given above, the spin 
current density must include two different parts in order to satisfy the generalized continuity equation. The first part 

is Jo + Jsoj it can be written as ^{s, jo -I- jso}- However, the jso cannot be simply written as jso = ^{ — ^^, (5(x — x)} 
(only the jo does) but should be written as a more complex expression given by Eq. psp . These two expressions are 



equal only if the velocity operator g' is a linear function of p (correspondingly, the degree of the full Haniiltonian 

must satisfies degiJ(p) < 2). The second part (Jlo)' is originated from the nonhnear spin-orbit coupHng terms, it 
appears only if deg/iso(p) > 2 and takes the form of commutator instead of anticommutator. From the viewpoint 
of noncommutability of operators, this new term is originated from the noncommutability of the density operator 
(5(x — x) and the ds/dt, because the latter one also includes momentum operators in spin-orbit coupling systems. 
This term does not give any additional contribution to the particle current (or charge current) but may contribute 
an additional surface integral term if we consider the time rate of change of the total spin in a given volume. 

As an example of these differences, we can consider a model with the following Dresselhaus spin-orbit coupling )^^'^^i^^ 
in which 

hso=l^<J^p^{pt+i-pf+2), {'i = x,y,z; i + 3^i), (44) 

i 

where a^ are the PauH matrices. After symmctrization, it can be rewritten as h^o — H!>JpiPjPk, where Hxyy — 

rr(3) _ A(3) _ _/.(3) _ _M3) _ _MS) _ 7^ 6^(3) _ £.(3) _ ^.(3) _ _ fri^) _ _ rT{3) _ _ rV(3) _ 7^ 
-iJ^yxy — ^^yyx — ^^xzz — ^^zxz — ^^zzx — 3 ^' V^^ — 2:yj3 — J-^zzy — ^^yxx — ^^xyx — ^^xxy — 3^y' 

Hfxx = Hxzx = Hxxz = -iifyy = ^ H yzy = ^ H^yz = ^(7 ^ , anA sll the other hI"',,^^ = 0. Because the Jo always 
has the same form for arbitrary model, so, here we only give the expression for 3 so, Jso ^^^ ^- According to the 
definitions (fTl), dUl) and Ull), we have 



(Jso). = l{{^x,{Dil^-&^J)p}-{<Ty,{D^^J+^'J)p} + {a.,{D^J+Dg^)p}) 



6 
6 

+PyPx + PxPy + PyPx)p} + {^z, {PxPz + Pzpx + Pxp'^ + PzPx + PxPz + PzPx)p})l' = l, (45) 



- {{^x, ify + PyPy + Py -p^- Pzp'z " Pz)p} - {^y , {PxPy + PyPx + PxPy 



(Jso)x = ^ {[ax, (Cf,) - C(/J)p] - [a,, (Cg) + (7(2))p] + [a,, {CfJ + C(/J)p] 

= g {[o-x, {Pl - py -Pl+ P'z)p] - [ay, {PxPy + PyPx + PxPy + P'yP'x)p] 

+ [0-z, iPxPz +PzPx+P'xP'z +PzPx)p\}l' = l, (46) 

The y and x components can be obtained by permutations of the suffixes according the rule x -^ y,y ^ z, z ^ x and 
a; — > z, — > y, 2/ — > a;, respectively. The spin torque density is 



^1^1 ^sS 



at, -^{P,+p'i) ((Pj+1 +P'i+if - (.Pi+2 + P',+2f) P 



(47) 



For this model, Jso can not be simply given by i{vso,<5(x — x)} with Vgo — dhso{p)/dp, although their results 
correspond to the same classical analogy. Because by the latter intuitive definition all terms such as p^p^, pxp'y 
etc. are replaced by p^Px (or p^p^), and pxPy (or p'^p'y), respectively, so it obviously violates the local equation of 

continuity. Moreover, the new term J^^, which is purely originated from the nonlinear spin-orbit interaction and 
cannot be given by intuitive method, must also be taken into account in order to satisfy the continuity equation. 

V. SUMMARY 

In conclusion, by defining a current and a torque density matrix, a general equation of continuity satisfied by 
the matrixes of the density, current density and torque density has been derived. This equation holds in spin-orbit 
coupling systems as long as their Hamiltonian can be expressed in terms of a power series in momentum. Thereby, the 
universal expressions of the current density matrix and torque density matrix have been uniquely determined. The 
current density matrix can completely describe the particle (or charge) current as well as the spin current in a unified 
form, while the torque density matrix represents spin precession caused by a total magnetic field, including a usual and 
an effective one. This effective magnetic field is originated from the spin-orbit coupling and depends on momentum. 
In contrast to the conventional intuitive definition, it is found that if the spin-orbit coupling Hamiltonian includes 
nonlinear term of momentum, then the definition of the current density matrix should contain a new additional 
term, which is purely originated from the nonlinear spin-orbit coupling and can only contribute to the spin current. 
Moreover, if the degree of the full Hamiltonian > 3, then the conventional intuitive definition of the current density 
must also be modified in order to satisfy the local conservation law of number. 
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APPENDIX A: DERIVATION OF EQ. (f20|) 

Because (1) pi is commute with f)[, and (2) i?j ..,; is symmetric under the interchange of its suffixes, we have 



H. 



j(m+l) 



-l\ rV-(™+l) p('") 



^":^:^. (pn • • -p.™.. - ^r^:^, J g(i, lo ^ ^ ^ - pd ^;::i^^.„ g(i, lo, 



(Al) 



where 



H. 



\{m+l) 



":i^. {pk ■ ■ -pL.,, - i€^:i,.) G(i, 1') == -2 fe -p:)if^™tiF':;:.:,,„G(i, lo, 



f,!^ A(m+1) p,(™) 



Fl 



- , - - d(1) I I d(™) 



(A2) 



(A3) 



F', 



(1) 



4---C+4---K™-i^»™+--- + < 



\(m) 



Therefore, we get 



H. 









where 



a 



(m) ^ p(m) _ p,(™) 



6(1) 









= Pil-- -Pirr. +Pil- ■ -Pi^-lR. 

Similarly, we have 

r. ^ /v\ ry(™+i) odC^+i) /»Ty(™+i) _ -*- /^ ^Mr'('") r'H-(™+i 

Adding Eq. ((J5|) and ((XT)) . we get 

By taking the limit as 1' ^ 1+, we obtain the Eq. (PO]) . 



/v(m+l) A(m) A 



ir{m+l) A(m+1) A 



(A4) 



(A5) 



(A6) 



(A7) 



(A8) 
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